Transport Phenomena and Structuring in Shear Flow of Suspensions near Solid Walls 
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In this paper we apply the lattice-Boltzmann method and an extension to particle suspensions 
as introduced by Ladd et al. to study transport phenomena and structuring effects of particles 
suspended in a fluid near sheared solid walls. We find that a particle free region arises near walls, 
which has a width depending on the shear rate and the particle concentration. The wall causes the 
formation of parallel particle layers at low concentrations, where the number of particles per layer 
decreases with increasing distance to the wall. 
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I. INTRODUCTION 

Many manufacturing processes involve the transport of 
solid particles suspended in a fluid in the form of slurries, 
colloids, polymers, or ceramics. Examples include the 
transport of solid material like grain or drug ingredients 
in water or other solvents through pipelines. Naturally, 
these systems occur in mud avalanches or the transport of 
soil in water streams. It is important for industrial appli- 
cations to obtain a detailed knowledge of those systems 
in order to optimize production processes or to prevent 
accidents. 

For industrial applications, systems with rigid bound- 
aries, e.g. a pipe wall, are of particular interest since 
structuring effects might occur in the solid fraction of 
the suspension. Such effects are known from dry gran- 
ular media restingon a plane surface or gliding down 
an inclined chute [l], Q . In addition, the wall causes a 
demixing of the solid and fluid components which might 
have an unwanted influence on the properties of the sus- 
pension. Near the wall one finds a thin lubrication layer 
which contains almost no particles and causes a so called 
"pseudo wall slip" . Due to this slip the suspension can 
be transported substantially faster and less energy is dis- 
sipated. 

The dynamics of single-particle motion, interaction 
with other particles, and effects on the bulk properties 
are well understood if the particle's intertia can be ne- 
glected. If massive particles are of concern, the behavior 
of the system is substantially harder to describe. A num- 
ber of people have studied particle suspensions near solid 
walls. These include Sukumaran and Seifert who describe 
the influence of shear flow on fluid vesicles near a wall , 
Raiskinmaki et al. who investigated non-spherical par- 
ticles in shear flow Jasberg et al. who researched 
hydrodynamical forces on particles near a solid wall 
and Qi and Luo who model the rotational and orienta- 
tional behaviour of spheroidal particles in Couette flows 
, or Ladd's work on the sedimentation of homogeneous 
suspensions of non-Brownian spheres . 

The last four authors use a simulation technique based 
on the lattice-Boltzmann equation (LBE), that we are 



also going to use in our simulations. 

Many other authors have studied similar systems the- 
oretically and experimentally. These include Chaoui and 
Feuillebois who performed theoretical and numerical in- 
vestigations on a single sphere in a shear flow close to 
a wall 0, IToL ITll ll^ l , or Datta and Shukla who pub- 
lished an asymptotic analysis for the effect of roughness 
on the motion of a sphere moving away from a wall [jj. 
Berlyand and Panchenko studied the effective shear vis- 
cosity of concentrated suspensions by a discrete network 
approximation technique 0,^3 an d Becker and McKin- 
ley analysed the stability of creeping plane Couette and 
Poiseuille flows [l5|. There is a theoretical and experi- 
mental study on rotational and translational motion of 
two close spheres in a fluid ^(|, and a general approach 
for the simulation of suspensions has been presented by 
Bossis and Brady 0, 0, 0] and applied by many au- 
thors. Melrose and Ball have performed detailed stud- 
ies of shear thickening colloids using Stokesian Dynamics 
simulations |20l l21| . Suspensions of asymmetric parti- 
cles like fibers, polymers, or large molecules have been 
of interest to many experimentalists and theoreticians, 
too. These include Schiek and Shaqfeh or Babcock et al. 

mm. 



We expect structuring close to a rigid wall at much 
smaller concentrations than in granular media because 
of long-range hydrodynamic interactions. In this paper, 
we study these effects by the means of particle volume 
concentrations versus distance to the wall. Autocorrela- 
tion functions of these profiles as well as autocorrelation 
functions of particle distances to a wall give detailed in- 
formation about the system's state and time dependent 
behavior. We study the dependence of correlation times 
on shear rates and achieve insight in the connection of 
the abovementioned lubrication layer on the shear rate 
and particle concentration. 

The remaining of this paper is organised as follows: 
After a description of the lattice-Boltzmann method and 
its extension to particle suspensions in the following sec- 
tion we give an overview about the simulation details in 
section ITTT1 Our results are presented in section Hvl and 
we conclude in section Ivl 
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II. SIMULATION METHOD 

The lattice-Boltzmann method is a simple scheme for 
simulating the dynamics of fluids. By incorporating solid 
particles into the model fluid and imposing the correct 
boundary condition at the solid/fluid interface, colloidal 
suspensions can be studied. Pioneering work on the de- 
velopment of this method has been done by Ladd et al. 
[24 125L |2a| and we use their approach to model sheared 
suspensions near solid walls. 



A. Simulation of the Fluid 

We use the lattice-Boltzmann (hereafter LB) simula- 
tion technique which is based on the well-established con- 
nection between the dynamics of a dilute gas and the 
Navier-Stokes equations [2^. We consider the time evo- 
lution of the one-particle velocity distribution function 
n(r, v, t), which defines the density of particles with ve- 
locity v around the space-time point (r, t) . By introduc- 
ing the assumption of molecular chaos, i.e. that succes- 
sive binary collisions in a dilute gas are uncorrelated, 
Boltzmann was able to derive the integro-differential 
equation for n named after him |27j 



dtn + v • Vn = 



(1) 



coll 



where the left hand side describes the change in n due to 
collisions. 

The LB technique arose from the realization that only 
a small set of discrete velocities is necessary to simulate 
the Navier-Stokes equations [28|. Much of the kinetic 
theory of dilute gases can be rewritten in a discretized 
version. The time evolution of the distribution functions 
n is described by a discrete analogue of the Boltzmann 
equation j2||: 

rn(r + CiAt,t + At) =m(r,t)+Ai(r,t), (2) 

where Aj is a multi-particle collision term. Here, rii(r,t) 
gives the density of particles with velocity c, at (r, t). 
In our simulations, we use 19 different discrete velocities 
Cj. The hydrodynamic fields, mass density p, momentum 
density j = pu, and momentum flux n, are moments of 
this velocity distribution: 

p=) j n i , j = pu = y^rtjCj, IT = rtjCjCj. (3) 



We use a linear collision operator, 

A i (r,t)=M ij (n j -n* q ) 



(4) 



where My = dA Q™. — - is the collision matrix and n eq 
the equilibrium distribution |29j, which determines the 



scattering rate between directions i and j. For mass and 
momentum conservation, M^- satisfies the constraints 



M 



M 



=0. 



(5) 



We further assume that the local particle distribution 
relaxes to an equilibrium state at a sing le rate t and 
obtain the lattice BGK collision term [3(j 



(6) 



A, = (m- n q ). 



By employing the Chapman-Enskog expansion |27], |3J| it 
can be shown that the equilibrium distribution 



nf = P uj c 



9 3 
l + 3c J -u+-(c l -u) 2 --u 2 



(7) 



with the coefficients of the three velocities 



1 



f 



,V5 



" - 3' 18' 
and the kinematic viscosity [2(| 

T) _ 2r - 1 

properly recovers the Navier-Stokes equations 



36' (8) 



du 
St 



1 7? 

(ttV)u = — Vp + - An, Vu = 0. 
P P 



Fluid-Particle interactions 



(9) 



(10) 



To simulate the hydrodynamic interactions between 
solid particles in suspensions, the lattice-Boltzmann 
model has to be modified to incorporate the boundary 
conditions imposed on the fluid by the solid particles. 
Stationary solid objects are introduced into the model by 
replacing the usual collision rules (Eq. ©) at a specified 
set of boundary nodes by the "link-bounce-back" colli- 
sion rule [i^]. When placed on the lattice, the boundary 
surface cuts some of the links between lattice nodes. The 
fluid particles moving along these links interact with the 
solid surface at boundary nodes placed halfway along the 
links. Thus, a discrete representation of the surface is 
obtained, which becomes more and more precise as the 
surface curvature gets smaller and which is exact for sur- 
faces parallel to lattice planes. Two discretized spherical 
surfaces near contact are shown as filled symbols in Fig. 
n Empty symbols denote the fluid, while filled squares 
and triangles depict the discretized surface. The crosses 
(C) denote the shared boundary nodes in contrast to the 
filled circle (E) which is not a shared boundary node since 
it is placed on individual links for each sphere. 
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Figure 1: The blank squares and triangles depict the inner 
fluid of particles A and B respectively. The filled squares and 
triangles denote the discretized surface; fluid particles cross- 
ing that surface are reflected. The empty circles represent 
the outer fluid, while crosses (C) denote the shared boundary 
nodes. The filled circle E is not a real shared node, but be- 
longs to both spheres (A and B), although it lies on different 
links. D is a single outer node between the two surfaces. If 
those surfaces move towards each other, high pressure occurs 
at D. 

Numerical results of simulations of a stationary 
Poisseuille-fiow between two flat surfaces are in good 
agreement with the theoretical formula |33| 

gL 2 ( 4(x-L) 2 \ , s 

This is demonstrated in figure [21 which shows the velocity 
profile vs. dimcnsionsless distance x/L from the left wall 
of a fluid with viscosity 77 = | and density p = 1 under 
constant force g exerted on each lattice point in a channel 
with width L. g is set to 10 -4 for L £ {8; 16; 32} and to 
g = 5 x 10~ 5 for L = 16. i],g,p,L are in lattice units. 
The solid line corresponds to the profile as expected from 
Eq. (ETJ. 

Since the velocities in the latticc-Boltzmann model are 
discrete, boundary conditions for moving suspended par- 
ticles cannot be implemented directly, fnstead, we can 
modify the density of returning particles in a way that 
the momentum transferred to the solid is the same as 
in the continuous velocity case. This is implemented by 
introducing an additional term Af, in Eq. (|2J) |24| : 

A M = 2ujC ' P f C \ (12) 

with c s being the velocity of sound and coefficients to Ci 
from Eq. JSJ . 

To avoid redistributing fluid mass from lattice nodes 
being covered or uncovered by solids, we allow interior 
fluid within closed surfaces. Its movement relaxes to 
the movement of the solid body on much shorter time 
scales than the characteristic hydrodynamic interaction 
[24l| . Fig. |21 shows a cut through a three-dimensional box 
containing a sphere S with periodic boundaries on front, 
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Figure 2: Poisseuille-flow of a fluid with viscosity 77 = | and 
density p = 1 under gravity g — 10 -4 exerted on each lattice 
point, channel width L £ {8; 16; 32}. Gravity is set to g = 
5 x 10 -5 for channel width L — 16. The solid line represents 
the expected profile (Eq. iflP ). 



back, left and right sides. On the top and bottom sides 
as well as on the sphere surface we use link-bounce-back 
boundary conditions. The particle is falling under the 
influence of gravity g. The system size is 32 x 32 x 32 
lattice constants a and the particle radius is 4a. At the 
beginning the particle and the fluid are at rest and after 
3000 time steps the particle attains a steady state. The 
cut in Fig. |3|has been generated after 5155 time steps, 
i.e. well after the system has reached the steady state. 
Its velocity u is 19% higher than u^, expected by Stokes' 
equation in an inifinite fluid system l33j | . The difference 
is caused by the fluid vortices V seen in Fig. |3 which is 
due to the periodic boundary conditions and could not 
arise in an infinite system. 



C. Boundary nodes shared between two particles 

If two particle surfaces approach each other within one 
lattice spacing, no fluid nodes are available between the 
solid surfaces (Fig. QJ. In this case, mass is not conserved 
anymore since boundary updates at each link produce a 
mass transfer At,a 3 (a =cell size) across the solid-fluid 
interface [24[ • The total mass transfer for any closed sur- 
face is zero, but if some links are cut by two surfaces, 
no solid-fluid interface is available anymore. Instead, the 
surface of each particle is not closed at the solid-solid con- 
tacts anymore and mass can be transferred in between 
suspended particles. Since fluid is constantly added or 
removed from the individual particles, they never reach 
a steady state. In such cases, the usual boundary-node 
update procedure is not sufficient and a symmetrical pro- 
cedure which takes account of both particles simultane- 
ously has to be used J2{|. Thus, the boundary-node ve- 
locity is taken to be the average of that computed from 
the velocities of each particle. Using this velocity, the 
fluid populations are updated (Eq. (|12f) ), and the force 
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Figure 3: Cut through a three-dimensional system after 5155 
time steps. The link-bounce-back boundary conditions are 
implemented on the surface of the sphere and the walls. The 
particle Reynolds number is Re = 6 and g is the gravitational 
force. The movement of the interior fluid has already relaxed 
to solid body movement. 

is computed; this force is then divided equally between 
the two particles. 

D. Lubrication interactions 

If two particles are in near contact, the fluid flow in 
the gap cannot be resolved by LB. For particle sizes used 
in our simulations (R < 5a), the lubrication breakdown 
in the calculation of the hydrodynamic interaction occurs 
at gaps less than 0.1R [32]. This effect "pushes" particles 
into each other. 

To avoid this force, which should only occur on in- 
termolecular distances, we use a lubrication correction 
method described in Ref . • For each pair of particles 
a force 

Fiub = -6tt?7 - - — U12 • | h<h K 

(i?i + R 2 y \h h N J |ri2 1 

(13) 

is calculated, where u 12 = u x — u 2 ,/i = |r 12 | — R\ — 
i?2 is the gap between the two surfaces and a cut off 
distance — fa 26]. For particle- wall contacts we 
apply the same formula with i? 2 — ► oo and h = |r 12 | — R\. 
The tangential lubrication can also be taken into account, 
but since it has a weaker logarithmic divergence and its 
breakdown does not lead to serious problems, we do not 
include it in our simulations. 

This divergent force can temporarily lead to high ve- 
locities, which destabilize the LB scheme. Instabilities 



can be reduced by avera ging the forces and torques over 
two successive time steps |25l|. In Ref. Q an implicit up- 
date of the particle velocity was proposed. This method 
then has been generalized and adopted for LB where two 
particles are in near contact |26l |3Jj ■ The drawback of 
this algorithm is the requirement of two sweeps over all 
boundary nodes. As we study creeping motion, we use 
the following simple method. High forces can only arise if 
the lubrication correction is switched on. Therefore, the 
lubrication correction Fj u b is limited to a value which 
would cause a particle acceleration of 0.1 Mach/s. Such 
a limitation may lead to particle overlap, but we found 
that on average there are only 5 occurences of this limi- 
tation per particle within 10 6 time steps. 

E. Particle motion 

The particle position and velocity is calculated using 
Newton's equations 

a = — F = v, v = r. (14) 

777 

The force F is obtained from the calculation of the 
particle-fluid coupling and the lubrication corrections. 
Then, the equations are discretized and integrated us- 
ing the Euler-Cromer method [35J. The velocity v n +i 
and position r„ + i for the time step n + l are obtained by 
utilizing the velocity, position and force from time step 
77 as well as the time step At = 1 and particle mass m. 

in 

v„+i = v„ + —At (15a) 

777 

r„+i = r„ + v„ +:L Ai (15b) 

The same method is applied to particle rotation, with 
position replaced by angles, velocity by angular velocity, 
force by torque and mass by moment of inertia. We do 
not use more sophisticated methods since they either re- 
quire additional memory and extra calculations (Verlet 
[36|. Runge-Kutta [33) or require the solution of an im- 
plicit equation for the velocity at each particle boundary 
node (Velocity- Verlet) (38|- Since the forces and veloci- 
ties in our simulation are rather small and the particle 
kinetic energy is not conserved between collisions (it is 
changed by particle-fluid interaction), we do not need to 
care for neglegible numerical inaccuracies of this method. 

III. SIMULATIONS 

The purpose of our simulations is the reproduction 
of rheologic experiment on computers. First, we simu- 
late a representative volume element of the experimental 
setup. Then we can compare our calculations with ex- 
perimentally accessible data, i.e. density profiles, time 
dependence of shear stress and shear rate. We also get 
experimentally inaccessible data from our simulations 
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like translational and rotational velocity distributions, 
particle-particle and particle- wall interaction frequencies. 

The experimental setup consists of a rheoscope with 
two spherical plates, which distance can be varied. The 
upper plate can be rotated either by exertion of a con- 
stant force or with a constant velocity, while the comple- 
mentary value is measured simultaneously. The mate- 
rial between the rheoscope plates consist of glass spheres 
suspended in a sugar-water solution. The radius of the 
spheres varies between 75 and 150 pm. For our simula- 
tions we assume an average particle radius of 112.5 pm. 
The density and viscosity of the sugar solution can also 
be changed. 

Because glass and suger solution have different light 
absorption constants, the particle concentration can be 
obtained by spectroscopic methods. Alternatively, the 
experimental material can be frozen and analyzed by an 
NMR spectroscope and a three dimensional porosity dis- 
tribution can be extracted from the data. Details of the 
experiment which is currently under development can be 

found in [nm mm. 

A low resolution (R ~ 2a) simulation of a system with 
the same volume as the experiment would need about 
10 GB RAM which is about five times as much as typ- 
ically available in current workstations. Each time step 
the program sweeps at least twice over the full data set. 
Simulating one minute real time would need about three 
years CPU time. Increasing the resolution or implement- 
ing curved boundary would increase the computation 
time even more. Therefore, we calculate only the be- 
havior of a representative volume element which has the 
experimental separation between walls, but a much lower 
extension in the other two dimensions than the experi- 
ment. In these directions we employ periodic boundary 
conditions for particles and for the fluid. 

Shearing is implemented using the "link-bounce-back" 
rule with an additional term Abj at the wall in the same 
way as already described for particles (Eq. i|12|) with 
now being the velocity of the wall). If a fluid node be- 
tween the particle and the wall is missing, we use the 
approach for shared boundary nodes as discussed in sec- 
tion EO 

To compare the numerical and experimental results, we 
need to find characteristic dimensionless quantities of the 
experiment which then determine the simulation param- 
eters. For this purpose we use the ratio of the rheoscope 
height and the particle size A, the particle Reynolds num- 
ber Re and the volume fraction of the particles <p: 

L T] V s 

with R being the particle radius, L the height of the rheo- 
scope, 7 the shear rate, r\ the fluid viscosity and N the 
number of particles. In the experiment the suspended 
particles have a slightly lower density than the fluid. Re- 
ducing the particle density would cause instabilities in 
LB. Therefore, we need to change the acceleration of 



gravity to a value, which would cause the same sedimen- 
tation or buoyancy velocity u. Stokes' law [3i| gives the 
connection between u and gravity g: 

F = 6ttRt]u u= — — ^— , (17) 

otvRt] 

with the effective mass m = ^TrR 3 (p s — p/) of the solid 
particle. Converting u to the dimensionless velocity u' 
(lattice constant/time step) and inserting simulation pa- 
rameters into the last equation we get 



where m' is the mass of the particle without interior fluid, 
R' the particle radius and rj the fluid viscosity (Eq. J5J)). 

To provide the simulation results with units, we cal- 
culate the length of the lattice constant a = R/R' 
and the duration of one time step At = 7'/7 : Using 
R = 1.125 ■ 10" 4 m, L = 3.375 ■ 10~ 3 m, p f = 1446 ^, 
p s = 1180 ■H, 17 = 0.450 7 = 10 s- 1 , R' = 

~° 5 ' m-s ' ' 30 ' 

L = 59, v' = | and p'± — 0.7 we obtain 

a = 0.572- 10~ 4 m, At = 1.262 ■ 10~ 4 s. (19) 

In the simulations presented in the next section we vary 
the particle Reynolds number to find the dependence of 
the time needed to attain a steady state and strength of 
structuring effects on the shear rate. 

Next we vary the particle volume fraction to study the 
correlation of velocity profiles and particle concentration. 
Different volume fractions lead to different correlation 
effects of particle positions and density profiles. 

To check our conversion rule between numerical and 
experimental data, we will try to change the fluid viscos- 
ity without changing the Reynolds number. This leads 
to different shear rates and consequently different time 
steps in the simulation. Higher viscosities lead to longer 
time steps and thus to shorter simulation times. 

A system with Re = 4-10 -6 needs about 900 s to attain 
the steady state. 900 s are equivalent to 7- 10 6 iterations. 
For such a high number of iterations the program requires 
about 20 CPU-days on a 2 GHz AMD Opteron. 

IV. RESULTS 

Fig. ^Ishows a snapshot of a suspension with 50 spheres 
after 5772500 time steps which are equivalent to 729 s. 
The vector g represents the direction of gravity and 
depicts the velocity of the sheared wall. 

The particles feel a gravitational acceleration g — 
0.8 m/s , have a mass m = 7.7T0~ 8 kg, a Reynolds num- 
ber Re = 4.066875- 10~ 4 , and a radius R = 1.125- 10" 4 m. 
The system size is 1.83 • 10~ 3 x 1.83 • 10~ 3 x 3.375 • 10~ 3 m 
which corresponds to a lattice size of 32 x 32 x 59. The 
density of the fluid is set to pf — 1446^ and its viscosity 
is rj = 450 mPa • s. The walls at the top and the bottom 
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m = 
2 and 



Figure 4: A snapshot of a suspension with 50 spheres (radius 
R = 1.125- 10" 4 m, mass m = 7.7- 10" 8 kg) at time t = 729 s 
which corresponds to 5772500 time steps. The volume of the 
simulated system is 1.83 ■ 10~ 3 x 1.83 • 10 -3 x 3.375 • 10 
11.3025 ■ 10~ 9 m 3 , acceleration of gravity g — 0.80 m/s 
shear velocity v s = 3.375 • 10 -2 m/s. The fluid has a viscosity 
r) = 450 mPa ■ s and density p/ = 1446 ^-§-. This visualiza- 
tion is a typical example for a system that has reached a 
steady state: All particles have fallen to the ground due to 
the exerted gravitational force and most of the system has no 
particles. 



are sheared with a relative velocity v s — 3.375- 10 -2 m/s. 
The system size, particle size and mass, as well as the 
gravitational force and all fluid parameters are fixed 
throughout the paper. After 200 time steps a linear fluid 
velocity profile can be observed and the particles are in- 
serted in a random fashion: After choosing a random po- 
sition for the particle, we check if the distance between 
these coordinates and the centers of all other spheres is 
at least 2R + a in order to avoid high interparticle forces. 
The initial particle velocities are set to the velocity of 
the fluid at the center of the particles. This algorithm 
allows a dense and uniform particle distribution within 
the whole simulation volume and has been applied in all 
our simulations. Fig. 0] is a representative visualization 
of our simulation data and demonstrates that after the 
system has reached its steady state, all particles have 
fallen to the ground due to the influence of the gravi- 
tational force. Most of the simulation volume is free of 



particles. 

In order to quantitatively characterize structuring ef- 
fects, we calculate the particle density profile of the sys- 
tem by dividing the whole system into layers parallel to 
the walls and calculating a partial volume Vij for each 
particle i crossing such a layer j. The scalar Vy is given 
by the volume fraction of particle i that is part of layer 
j- 



R 2 (R 



max 
ij 



nnun 
H ij ) 



max 



RT)) (20) 



If the component z perpendicular to the wall of the 
radius vector of the center of sphere i lies between 
we have 



r mm and r n 



=[j--)AL z -R, 



= [j + -)AL z + R, 



and 



Din ax 

H ij - ' 


[R 

max . 


if n iZ 
else 


+ R < rf 3 * 


nmin 


l-R 

„min 


if n,« 
else 


— R> rf n 



Finally, the sum of all weights associated with a layer is 
divided by the volume of the layer 



1 



N 



M 



(21) 



with L x , L y being the system dimensions between peri- 
odic boundaries, L z the distance between walls, M the 
number of layers, and AL Z the width of a single layer. 

Density profiles calculated by this means for systems 
with two different shear rates 7 = 10 s _1 and 7 = 1 s _1 
are presented in Fig. All other parameters are equal 
to the set given in the last paragraph. The peaks in Fig. 
^demonstrate that at certain distances from the wall the 
number of particles is substantially higher than at other 
positions. The first peak in both figures is slightly below 
one particle diameter, which can be explained by a a lu- 
bricating fluid film between the first layer and the wall 
which is slightly thinner than one particle radius. Due 
to the small amount of particles, time dependent fluc- 
tuations of the width of the lubricating layer cannot be 
neglected and a calculation of the exact value is not pos- 
sible. The five peaks in Fig. [SJi have samilar distances 
which are equal to one particle diameter. These peaks 
can be explained by closely packed parallel layers of par- 
ticles. Due to the linear velocity profile in z direction of 
the fluid flow, every layer adopts the local velocity of the 
fluid resulting in a relative velocity difference between 
two layers of about 2i?7. These layers stay stable in time 
with only a small number of particles being able to be 
exchanged between them. 
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Figure 5: Density profiles from simulations with two different shear rates 7 = 10 s" 1 (a) and 7=1 s" 1 (b). Other parameters 
are equal to those given in Fig. 0] (a) shows five peaks with separations about one particle diameter, which reveal the forming 
of particle layers. The number of particles per layer is decreasing with increasing distance to the wall, and the change in particle 
numbers is caused by gravity which is directed perpendicular to the wall at z = 0. Although we used the same gravity and 
particle numbers, there are only three peaks in (b) and their width is higher than in (a), demonstrating that the structuring 
effects strongly relate to the shear rate. 



Fig. |3> only shows three peaks with larger distances 
than in Fig. [JJl However, the average slope of the pro- 
file is identical for both shear rates. For smaller shear 
rates, velocity differences between individual layers are 
smaller, too. As a result, particles feel less resistance 
while moving from one layer to another. Every inter- 
layer transition destorts the well defined peak structure 
of the density distribution resulting in only three clearly 
visible peaks in Fig. [SJd. 

With changing time, the first peak stays constant for 
both shear rates. The shape, number and position of all 
other peaks is slightly changing in time. 

To aquire a quantitative description of ths effect we 
calculate the autocorrelation function of the density pro- 
file (Fig. H3) r T for each individual layer I, 

(rqjE 0'(j--Ai)-^((i + j)At) 
r l T (i-M) = ^ = (22) 

j'=i 

with At being the time step, i the current iteration and 
T the total number of time steps. 

Averaging the r\ over all M layers gives 

1 M 

r T (i-At) = -^2r l T (i-At), (23) 
1=1 

which is presented in Fig. 0for two systems with shear 
rates 7 = 10 s _1 and 7 = 1 s _1 . The autocorrela- 
tion starts — as given by definition — at one. Then, 
it decreases and converges to constant values at about 
i ■ At = 15 s for 7 = 10 s _1 and at i ■ At = 25 s for 
7 = 1 s" 1 . We obtain these values by fitting the data 
to a constant function using a nonlinear least-squares 
Marquardt-Levenberg algorithm. The computed values 



of the autocorrelation function are different for the given 
shear rates: r T is 0.480 for 7 = 1 s _1 and 0.361 for 
7 = 10 s^ 1 respectively. It is evident that for a simula- 
tion without shear the autocorrelation converges to one 
because after sedimentation the density profile should 
not change. Thus, (j> l (k ■ At) is almost constant for all 
fc, and r> — > 00. For 7 — > 00 the velocity and the col- 
lision frequency are increasing and the correlation de- 
creases for high shear rates. Therefore, the expectation 
that for smaller shear rates the autocorrelation converges 
to higher values than for larger shear rates, is confirmed. 

Another possibility to compute typical correlation 
times of structured layers is to analyze the autocorre- 
lation of particle distances to one of the walls. For this 
purpose we replace the volume fraction 4> l (k ■ At) of layer 
I by the distance of particle / to one of the walls r z in 
Eq. (|22|l . Then the acquired data is averaged for all N 
particles 

1 N 

r r (i.At) = -£V T (i.Af). (24) 
1=1 

The dependence of r T on time calculated by this means 
is shown in Fig. [7| where simulation parameters are as 
in the previous section. It is possible to fit the data to 
an exponential function of the form 

r r=e -^7, (25) 

where r corr is the characteristic correlation time. We 
get T corr = 5.5 s and r corr = 38.64 s for 7 = 10 s _1 
and 7 = Is -1 , respectively. This fully corresponds to 
the behavior expected from the density profiles: Shorter 
correlation times are related to higher shear rates. At 
higher shear rates the mean velocity of the particles is 
also higher. Thus they collide with other particles and 
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Figure 6: Autocorrelation r T of the density profiles shown in Fig. In both plots the autocorrelation converges to a fixed 
value. The dashed lines correspond to fitted constants and the points to the simulation data. In (a) the system has a shear rate 
7 = 10 s _1 and in (b) 7 = 1 s _1 . The higher shear rate leads to higher particle velocities and therefore to a higher collision 
frequency. Therefore, this system faster attains the steady state and is less correlated. This is confirmed by the smaller limit 
of r T , which is 0.361 instead of 0.480 in (b). 
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Figure 7: Autocorrelation r T of the particle distances from a wall for two systems with shear rates 7 = 10 s 1 (a) and 7 = 1 s 1 
(b). All simulation parameters exept 7 correspond to those given in Fig. 01 The straight line in the plots with logarithmically 

The typical correlation time r evidently 



scaled r T -axes shows the exponential connection of r T and time t: r T oc e *^ Tc ' 
depends on the shear rate. We find r cor r ~ 5.5 s and r corr = 38.64 s for 7 = 10 s~ 



and 7 = 1 s respectively. 



walls more often. Each collision contributes a random with 
uncorrelated force component to the equation of motion, 
which reduces the correlation of particle positions. 

We also expect a strong dependence on the average 
particle concentration and different values for the grav- 
itational force. For a larger number of particles in the 
system, the effective viscosity changes which influences 
the collision rate and reduces correlation times. 

Also for very high shear rates there should be nonzero 
correlation times, and we expect a non-linear connection 
between shear rate and correlation time. Therefore, we 
did the same calculations for more different shear rates 
and plot the correlation times r corr versus shear rates 7 
in Fig. [SJ Rescaling the axis of ordinates logarithmically 
we obtain a straight line again. For high shear rates the with 
correlation time is decreasing exponentially: 

W = "C?-e-£, (26) 



: = 47.24 s being the maximum correlation time 
and 70 = 4.78 s _1 being a characteristic shear rate. 

Another interesting property is the distribution of par- 
ticle distances, which can be acquired by calculating the 
distances r(i,j, I) of all particle pairs. Because of periodic 
boundary conditions we also account for particle pairs if 
one of them is shifted in one of the nine possible periodic 
directions. The maximum distance is then limited by the 
smallest system dimension L m [ n (here 1.83 • 10 -3 m). 



, N N 9 (- 1 

p k = —Y y y if 

w ' LLZj '0 else 



r(i,j,l) 
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(27) 
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Figure 8: Dependency of the correlation time r corr on the 
shear rate 7. All data points lie on a straight line, which 
indicates an exponential behaviour of r CO rr: T CO rr oc e~ 
with 70 = 4.78 s" 1 . 
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Each pk corresponds to an inter-particle distance 



k+ - ) Ar. 



(28) 



In a homogeneous system the number of particles with a 
distance 6 to a given particle is proportional to the sur- 
face of a sphere of radius 5. Thus to avoid overweighting 
of larger particle distances we divide the number of par- 
ticles with distance 6 by 6 2 (Eq. (|2*7jl ). 

In Fig. ®l we present two distributions for a system 
with 50 particles. The first distribution corresponds to 
the start of the simulation at t = s. It has one peak 
between 6 = 1 and 6 = 2, after which it decreases con- 
tinuously. The measurement in steady state gives the 
second distribution at t = 865 s. This distribution also 
has one peak, but it is narrower and much higher than at 
t = s. The position of this peak corresponds to a dis- 
tance 6 slightly higher than 1 particle diameter to each 
other, i.e. most particles have a distance of about one 
particle diameter. 

By computing only the distribution of the components 
of particle distances perpendicular to the walls r z for the 
same system we get the results plotted in Fig.l^Ja. We do 



not need to account for periodic boundaries here resulting 
in a smaller number of counted particle pairs and slightly 
worse statistics: 



jv 



n 



Pk 



N(N 



' i=i j=i+i 




(29) 



with N being the number of particles in the system, and 



A? 



M ' 



(30) 



where M is the resolution of the distribution, L z the dis- 
tance between the walls, and 6 is calculated as given in 
Eq. J3SJ). The distributions in Fig. 03 show that for 
t = s there are no structured layers. This histogram 
gives a nearly straight line with a negative slope. 

Let us consider a homogeneous system completely 
filled with spheres and x being the number of particles 
in an individual layer. Then, x(x — 1) is the number of 
particle pairs with distance 5 = L z — 1. For 5 = L z — 2 
there are two pairs of layers of that distance. Thus, we 
get 2x(x— 1) particle pairs. Reducing the distance by one 
particle diameter increases the number of possible parti- 
cle pairs by x(x — 1). The total number of particle pairs 
is propotional to L z — 6. This argumentation is valid 
for all homogenously filled systems. This consideration 
is comfirmed by Fig. ^Jjp for t = 0. The slope of the line 
should be propotional to the volume fraction because x 
gets larger for higher particle numbers. After 865 s there 
is a peak for 6 = showing that most pairs belong to the 
same layer. The second peak belongs to 6 ~ 1, i.e. to 
particle pairs at adjacent layers. 

Both histograms in Fig. have a clear dependence 
on time. To visualize that dependence we calculate the 
mean (6) and standard deviation a of particle distances. 



1 M 



\ 



1 M 

-J2p„-(S 1 ~(6))\ 

i=l 



(31a) 
(31b) 



6, = 



i ■ 6 n 



M 1 



with M being the resolution of the distribution p(S), and 
<5max being the maximum particle distance. For the dis- 
tribution of differences of particle distances to one of the 
walls <5 max is set to 3.375T0 -3 m = 15 particle diameters. 
For the distribution of particle distances <5 max = 1.8 ■ 
10~ 3 m = 8.136 particle diameters. In Fig. HOfa the 
mean (6) (left ordinate) and a (right ordinate) are plotted 
over time. For short times, the mean decreases almost 
linearly to t ~ 16 s, then the slope approaches and (6) 
fluctuates around 6 = 4.65 till the end of the simulation 
(t ~ 900 s). At the beginning of the simulation it is not 
possible to recognise the characteristics of the evolution 
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Figure 9: Distributions of particle distances (a) and of differences of particle distances to one of the walls (b). All simulation 
parameters are equal to those given in Fig. [I] In both figures we show histograms for two different times: t — s and t — 865 s. 
In (a) the dominant peak moves from 5 ~ 1.2 to S ~ 1 showing the compression of the system under the influence of gravity. 
There is also a shoulder on the right of the peak showing that many particles have distances between 1 and 2 particle radii. In 
(b) we see a linear profile for t = 0, caused by the homogeneous particle distribution at the beginning of the simulation. The 
highest peak is at 8 — 0, which is caused by particles belonging to identical layers. The following peak is due to particles from 
adjacent layers. 
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Figure 10: The mean (6) and standard deviation a of particle distances (a) and differences of particle distances to one of the 
walls (b) versus time for a system of 50 particles with radius R = 1.125 • 10 -4 m, acceleration of gravity g = 0.80 m/s 2 . The 
Reynolds number is Re = 4.066875 • 10 -4 , and the shear rate 7 = 10 s _1 . In both figures (<$) and a converge to specific values 
within ca. 15 s and they fluctuate around fixed values till the end of the simulation (t ~ 900 s). In both cases the standard 
deviation is smaller than the mean. The mean particle distance converges to (6) ~ 4.6 and its standard deviation converges to 
a ~ 3.2. The mean difference of particle distances to the wall converges to (a) ~ 1.2 and a ~ 1. 



of a. For times between 6 and 15 s the points of a lie 
nearly on a straight line. At t ~ 16 the slope of a be- 
comes zero and a is fluctuating around a mean a = 3.2 
like (S). 

In Fig. 1 10b the mean (S) and standard deviation a 
of differences of particle distances to one of the walls 
are plotted versus time. To calculate these values we 
used the equations ("31J1 . The evolution of (S) and a is 
nearly linear between t — 5 s and t = 12 s. The slope 
then vanishes and only some random fluctuations can be 
seen around (6) — 1.4 and a — 1.2 particle diameters for 
t > 17 s. Note that the particle distances attain a steady 
state already after 15 s while the density profile needs 
158 s. 



To study the demixing phenomena already demon- 
strated in Fig. [SI we analyze the dependency of the 
particle and fluid velocities on the distance to the wall. 
Both profiles in Fig. II II are for a system with shear rate 
7=10 s _1 at t = 865 s. All other simulation parameters 
are kept as in the last section. In addition to the velocity 
profiles we plot a solid line corresponding to the fluid ve- 
locity profile of a system without particles. The values of 
fluid velocities at the walls (z = and z — 15 particle di- 
ameters) exactly match the wall velocities: v(0) = and 
u(15) = 3.375 • 10~ 2 m/s. For 2 < z < 6 both profiles 
agree very well with each other. No particles are present 
above z = 6 and the fluid velocity profile is exactly lin- 
ear. We do not have any particle velocity data in this 
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(a) z [particle diameter] (b) z [particle diameter] 

Figure 11: Velocity profile of a system with shear rate 7 = 10 s _1 and mean volume fraction (f> — 0.026 versus distance to 
one of the walls z. Solid lines correspond to the expected fluid velocity profile in absence of particles. At the walls (2 = 
and z = 15), the fluid velocity is identical to the wall velocities and in the particle filled region the fluid and solid velocities 
are equal confirming the validity of the no-slip boundary conditions on particle and wall surfaces. The velocity of the solid 
particles does not disappear at the wall unlike the fluid velocity, but converges to a fixed value instead. 



case. Below z — 2 the profiles separate and for z < 0.5 
the fluid velocity profile corresponds to the expected pro- 
file for a particle free system, while the particle velocities 
stay constant. This can be seen in Fig. Illb . which shows 
an enlarged particle velocity profile. The particle veloc- 
ity converges to v s = 1.1 • 10~ 3 m/s for z — > 0. For 
higher values of z it is linear, but its slope is about 10% 
smaller than the slope of the solid line. For z > 6 the ve- 
locity profile is linear, but it raises faster than expected 
in order to fit the wall velocity at z = 15 and to con- 
serve the validity of the no-slip boundary conditions at 
the walls. Since the particle and fluid velocities are iden- 
tical for 2 < z < 6, the no-slip boundary conditions on 
the particle surface are shown to generate correct results, 
too. 

The dependence of the particle velocity near the wall 
on the shear rate is studied in Fig. Illb for 7 = 1,0.1, and 
0.25 s _1 by calculating the particle velocities for z — » 0. 
Fig. ^Jdepicts these velocities versus shear rate and their 
linear dependence is clearly observable. The data from 
simulations with higher particle concentration (4> = 0.053 
instead of <j> = 0.026) also gives a straight line but with 
smaller slope. The slopes can be interpreted as an effec- 
tive width of the particle free layer near the wall, which is 
is 1.16-10" 4 m or 9.23-10" 5 m for <p = 0.053 or <j) = 0.026, 
respectively. The value for <fi = 0.053 is slightly smaller 
than the particle radius and in good agreement with ob- 
servations from the particle concentration profiles in Fig. 
[5^. The smaller width of the particle free layer at higher 
particle concentrations is caused by the higher pressure 
on the lowest particle layer. Since the system cross sec- 
tion is the same in both simulation series, with higher 
particle number the number of the particle layers in- 
creases. Thus, the resulting gravitational force on the 
lowest layer increases proportionally to the particle num- 
ber. However, the reciprocal width of the particle free 
layer is not proportional to the particle number because 
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Figure 12: The velocity of the "pseudo- wall-slip" versus shear 
rate for two different volume concentrations. The dependence 
of the velocity is linear. The slope of the line gives an effective 
width of the particle free region near the wall. The width 
is 1.16 ■ 10" 4 m and 9.23 • 10 -5 m for cf> = 0.026 and <j> = 
0.053, respectively. Narrower particle free regions are caused 
by higher forces due to weight of particles being above the 
particle layer near the wall. 



this layer is caused by the competition of gravity and the 
resistance to particle motion perpendicular to the wall. 
This is not constant but rather approximately propor- 
tional to the reciprocal value of the distance |4i| • 

We calculate the distributions of velocity components 
in three directions: Perpendicular to the wall (Fig. 113b). 
parallel to the shear direction (Fig. ^] perpendicular to 
the shear direction and parallel to the wall (Fig. 113b) 
In Fig. E| one can clearly recognise three peaks. The 
first peak is at 1.1 • 10~ 3 m/s, exactly corresponding to 
the wall slip velocity for the given shear rate. It can 
be seen that this peak corresponds to the lowest particle 
layer and all three peaks have a distance of 2.5 • 10~ 3 m/s 
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Figure 13: Distributions of particle velocities averaged over 5.55 • 10 time steps of the steady state, (a) shows the component 
perpendicular to the wall, and (b) perpendicular to the shear velocity. While the means of both velocity distributions are zero, 
their widths differ. The movement perpendicular to the wall is restricted by walls and structured layers. Both data cannot be 
fitted by a Gaussian distribution. 
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smaller particle numbers per layer their movement within 
the layer is less restricted resulting in the possibility to 
achieve higher intcr-laycr particle velocities. 

Particle velocity distributions perpendicular to the 
wall and parallel to the wall but perpendicular to the 
shear direction are presented in Fig. 113b and 113b re- 
spectively. The means of both distributions are zero as 
expected. The distribution of particle velocities perpen- 
dicular to the wall is narrower because the movement 
to the wall is restricted by lubrication interactions. The 
change between the layers is restricted by the differences 
in layer velocities, but it is not completely impossible. 
The data of both distributions do not follow a Gauss- 
distribution. 



Figure 14: Distribution of the particle velocitiy compo- 
nent parallel to the shear direction averaged over 5.55 • 10 6 
time steps of the steady state. The peaks are separated by 
2.5 • 10 -3 m/s, starting at 1.1 • 10~ 3 m/s. Dividing this veloc- 
ity difference (2.5 ■ 10 -3 m/s) by the shear rate (i.e. 10 s -1 ) 
results in the particle diameter since the average width of the 
layer corresponds to one particle diameter. These layers move 
against each other with a relative velocity corresponding to 
the shear rate. 



which matches the product of the shear rate and particle 
diameter. We have already seen the formation of particle 
layers near the wall, with a distance of about one parti- 
cle diameter (Fig. [SJ). Therefore, we assume that each 
peak in Fig. ^] is caused by one single particle layer. 
The height of the peaks decreases with the velocity since 
the number of particles per layer is being reduced with 
time (see Fig. |SJ. This reduces the probability of finding 
a particle with the velocity of the layer, which on the 
other hand is decreasing with the distance to the wall 
(Fig. Ill|) . Thus, for higher wall distances we get higher 
velocities and smaller particle numbers. The width of the 
peaks in Fig. ^] is increasing with the velocity. Due to 



V. CONCLUSION 

We successfully applied the lattice Boltzmann method 
and its extension to particle suspensions to simulate 
transport phenomena and structuring effects under shear 
near solid walls. We adopted the simulation parameters 
to the experimental setup of Buggisch et al. [33 and are 
able to obtain not only qualitatively comparable results, 
but also values that quantitively correspond to experi- 
mentally measured parameters. We hope to be able to 
report on direct comparisons between our theoretical re- 
sults and the experimental results of Buggisch et al. in 
the near future. 

We have shown that the density profile has several 
peaks, confirming the formation of particle layers. The 
density profile is changing in time, but its autocorrela- 
tion function converges to a non-zero value. On the other 
hand the autocorrelation function of particle distances to 
a wall converges exponentially to zero resulting in a fixed 
correlation time. This time is exponentially depending 
on the shear rate. Furthermore, we have shown that the 
particle distances attain a steady state at a much ear- 
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lier state of the simulation than the density profile. We 
have also shown the occurrence of a "pseudo- wall-slip" of 
particles, exhibited by a particle free fluid layer near the 
wall. The velocity of this slip has a linear dependence 
on the shear rate. It is possible to calculate an effective 
width of the particle free layer, which depends on the 
particle concentration. 

A natural extension of this work would be to increase 
the size of the simulated system in order to reach the 
dimensions of the experimental setup. Even though the 
number of LB time steps of our simulations is extremely 
high already, even longer runs would be desirable. 

It would also be interesting to study the behavior of 
the system for higher particle densities and higher shear 
rates. However, improvements of the method are manda- 
tory in order to prevent instabilities of the simulation. 
Without further improvement of the simulation method, 
the maximum particle volume concentration is limited 
to 0.3 and the maximum available shear rate is about 



10 s _1 . A possible solution of this well-known problem is 
the replacement of the velocity update by an implicitc 
scheme [2^, [3^. The artifacts caused by the interior 
fluid can be removed by slightly modifying the coupling 
rules 0| . We have not implemented this because of the 
high numerical effort, which is caused by the necessity to 
sweep over all boundary nodes twice, in order to redis- 
tribute the mass from nodes being covered by the spheres. 
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